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ABSTRACT 

N-body simulations suggest that the substructures that survive inside dark matter 
haloes follow universal distributions in mass and radial number density. We demon¬ 
strate that a simple analytical model can explain these subhalo distributions as re¬ 
sulting from tidal stripping which increasingly reduces the mass of subhaloes with 
decreasing halo-centric distance. As a starting point, the spatial distribution of sub¬ 
haloes of any given infall mass is shown to be largely indistinguishable from the overall 
mass distribution of the host halo. Using a physically motivated statistical description 
of the amount of mass stripped from individual subhaloes, the model fully describes 
the joint distribution of subhaloes in final mass, infall mass and radius. As a result, 
it can be used to predict several derived distributions involving combinations of these 
quantities including, but not limited to, the universal subhalo mass function, the sub¬ 
halo spatial distribution, the gravitational lensing profile, the dark matter annihilation 
radiation profile and boost factor. This model clarifies a common confusion when com¬ 
paring the spatial distributions of galaxies and subhaloes, the so called “anti-bias”, as 
a simple selection effect. We provide a Python code SubGen for populating haloes 
with subhaloes at http; //icc. dur. ac. uk/data/ 
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1 INTRODUCTION 

With advances in numerical simulations, the statistics of 
subhaloes inside dark matter haloes are being quantified 
to ever higher accuracy. These statistics cover a variety of 
different properties of subhaloes, including mass, position, 
size, maximum circular velocity, orbit, mass stripping rate 
and accretion time. Among these statistics, two outstand¬ 
ingly simple distributions are the subhalo mass function and 
the spatial distribution. It is well established that the mass 
function of subhaloes follows a universal power-law form 
(except for an exponential high mass tail), dA/dlnm = 
AMhostm-““, with a universal amplitude. A, and a univer¬ 
sal slope, a ~ 0.9, t hat are both independent of the host 
halo mass, Mhost (e.g. lOao et al]l2004bl . l2012bl : lOiocoli et al.l 
l2008h . The spatial distribution is known to be less con- 
cent rated than both t he DM and galaxy d i stributions (see 
e.g., Gao et aDl2004al : Ibibeskind et al.ir2005l : IPiemand et al.l 
I 2 OO 4 . and references therein). With recent zoom simula¬ 
tions that are able to resolve subhaloes spanning ~ 7 orders 
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of magnitude in mass and ~ 2 orders in separation from 
the centre of the host halo, the spatial distribution of sub¬ 
haloes is also found to be universal, that is independent of 
the mass of subhalo l|Springel et al.l[2008al : [^ao et al.ll2012bl : 

iHellwing et ahllioisl). 


Tracking progress in simulations are models for 


of the models are semi-analvtical (e.g.. 

Tavlor & Babull 

2 OOII: iBenson et al.l 

2 OO 2 I: iTavlor & Babul 

|2004 2005albl; 

Zentner et al.l 20051: 

Penarrubia & BensonI 

2 OO 5 I). typically 


starting from Monte-Carlo merger trees to evolve subhaloes 
dynamically after infall. Each subhalo is assigned a set of ini¬ 
tial orbital parameters according to distributions extracted 
from simulations. The orbits of these subhaloes are then 
evolved inside the host halo, with the mass of the subhalo 
updated at each timestep according to the tidal radius and 
a dynamical timescale. The statistical properties of the final 
subhalo population, such as the mass and velocity functions 
and the spatial distribution can be obtained. These models 
typically involve a number of free parameters that are cali¬ 
brated by comparing the predicted subhalo distributions to 
simulations. The advantages of such models are that they 
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can incorporate detailed physical processes, including dy¬ 
namical friction, tidal heating, tidal stripping and total dis¬ 
ruption. Furthermore, these models are often incorporated 
in semi-analytic models of galaxy formation. However, being 
a semi-analytical model means one has to resort to Monte- 
Carlo realizations of merger trees and detailed evolution of 
individual orbits to obtain population statistics. 

_ A simplihcation was introduced bv ivan den Bosch et al.l 

ll2005h . in which they only considered the average mass loss 
rate of subhaloes without evolving individual orbits. With 
information on the infall mass and infall time of subhaloes 
obtained from Monte-Carlo merger trees, they can evolve the 
overall mass distribution of subhaloes to recover the final 
mass and velocity functions (I Jiang fc van den Boschll2014l : 
Ivan den Bosch fc Jiang|[2014l l. However, the price of not fol¬ 
lowing individual orbits is the loss of spatial information, so 
that these models cannot predict the spatial distribution of 
subhaloes. 

In this work, we propose a simple and fully analytical 
model that simultaneously predicts both the mass function 
and spatial distribution of subhaloes, focusing on the key 
ingredients that shape these distributions. We start from 
empirical results on the infall mass function and spatial dis¬ 
tribution of subhaloes labelled by their infall masses (i.e., 
mass at accretion), which are found to be simple to describe 
and easy to interpret. The only difference between these dis¬ 
tributions and the final distributions is reduction in mass of 
each subhalo. For this reason, we will call the distribution 
labelled by infall mass as “unevolved” and that labelled by 
hnal mass as “evolved”. To obtain the latter, we only need 
to specify the connection between the final mass and the 
infall mass of each subhalo. This is achieved by a physically 
motivated statistical description of final-to-infall mass ratio 
at each halo-centric radius, rather than by relying on recipes 
to evolve the mass and orbit of each subhalo individually. 

The idea that the mass and spatial distribution can 
be derived from coupling the unevolved di stributions with 
subhalo stripping w as already explored by iLed ll2004l ') and 
lOguri fc Led (I 2 OO 4 II. An extended Press-Schechter (hereafter 
EPS; jBardeen_jtajJ jl986l : iBond et"^ Il99ll : iLacev fc Cold 
Il993l : l^ieth et al.l2001^ approach was adopted to predict the 
progenitor mass function at each radius inside the host. The 
mass and location of these progenitors were then adjusted 
following simple tidal stripping and dynamical friction pre¬ 
dictions. The resulting mass and spatial distributions of sub¬ 
haloes largely agreed with the results from N-body simula¬ 
tion available at that time. These models are fully analytical, 
with only one free parameter in iLed (l2004h a nd essentially 
zero free parameters in the improved version (lOguri fc Led 
|20041 . However, the model assumptions are not directly val¬ 
idated by simulations, and some of them are too idealized to 
be realistic. For example, subhaloes are all placed on circular 
orbits, which is certainly not the case in cosmological simu¬ 
lations. They also assume that initially the subhalo spatial 
distribution traces the density profile of the host halo. This 
differs from the unevolved spatial distribution (the distribu¬ 
tion of subhaloes selected by their infall mass) in that the 
former describes the location of the progenitor haloes at the 
formation time, while the latter describes the location at the 
hnal time. However, the circular orbit assumption, coupled 
with the assumed initial spatial distribution, leads to the 
same unevolved spatial distribution as we hnd in this work. 


This could explain why they are largely able to reproduce 
the subhalo spatial distribution in simulations despite unre¬ 
alistic assumptions. In contrast to their models which mostly 
start from theoretical assumptions, our model is built upon 
empirically validated assumptions. The success of our model 
thus serves as a promising starting point for more realistic 
hrst-principle models that attempt to understand fully the 
formation and evolution of subhaloes. 

We use two sets of hi gh-resolution zoom-in simula¬ 
tions (the Aquari us project, S pringel et al.|[2008al . and the 
Phoenix project, JCao et al.ll2012aH to verify and calibrate 
our model. The high resolution of these simulations enables 
us to make a detailed statistical analysis of the spatial distri¬ 
bution and stripping of subhaloes with unprecedented pre¬ 
cision. Each halo is also simulated at a series of resolutions, 
allowing us to distinguish physical properties from numerical 
artefacts by carrying out convergence studies. 

With the calibrated assumptions, the model simulta¬ 
neously recovers both the hnal subhalo mass function and 
spatial distribution accurately. The full model specihes an¬ 
alytically the joint distribution of subhalo infall mass, hnal 
mass and location inside the host. Such a joint distribution 
is ideal to obtain fast realizations and for Halo Occupation 
Distribution (HOD) modelling of subhaloes inside a host 
halo. We provide a Python implementation, SubGen, for 
such purposes. We also present several example applications, 
including the implication to the universality of the subhalo 
mass function, a prediction for lensing measurements of sub¬ 
halo masses and the modelling of dark matter annihilation 
radiation including the boost factor from subhaloes. 

This paper is organized as follows: the simulations used 
to verify and calibrate the model are introduced in Section[2] 
the key idea and the ehectiveness of our model is briehy 
demonstrated using a simplihed version in Section O the 
three model assumptions are then investigated and validated 
in detail in Section [d] the full statistical model is presented 
and extended to arbitrary host masses in Section [S] some 
implications and applications of the model are discussed in 
Section |6l hnally we summarize and conclude in Section [3 

In most cases, we use lower case to refer to properties of 
the subhaloes, and capitals to refer to the properties of the 
host. For example, m and r refer to the mass and radius of 
a subhalo, and M and R refer to the mass of the host and 
the radial location within a host. The virial mass, M200, is 
defined to be the mass inside the virial radius, R200, inside 
which the average density equals 200 times the critical den¬ 
sity of the universe. We use a mass scale mo = 10^°h“^MQ 
to normalize the halo mass in scaling relations. 


2 SIMULATIONS 

We use two sets of high resolution simulations to ver- 
ify and calibrate our model assumptions. One is Aquarius 
llSpringel et al.l^2008a^ . a set of zoomed-in simulations of 
Milky Way sized haloes each with a mass of ~ lO^^/i“^M 0 . 
Six different haloes are simulated (named A to F) in total, 
at hve different levels of resolutions (labelled 1 to 5 from 
high to low resolution, with particle masses ranging from 
~ lO^h~^ M 0 to ~ lO®/i“^M 0 ). The other is the Phoenix 
llGao et al.ll2012fl l simulations of nine cluster sized haloes 
each with a mass of ~ lO^®/i“^M 0 , run at four different 
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levels of resolutions (labelled 1 to 4, with particle masses 
ranging from 6 x to so that the same 

level number corresponds to a similar number of particles 
inside the host halo both in the Aquarius and Phoenix sim¬ 
ulations). For demonstartion purposes we will mainly focus 
on Aquarius halo A (named AqAl to AqA5 from level 1 to 
5). The properties studied in this work behave qualitatively 
similar across the different haloes, with only small differ¬ 
ences in parameter values. Both simulations are run with 
cosmological parameters = 0.25, Ha = 0.75, as ~ 0.9, 
ris = 1 and h — 0.73, where Hq = lOOh km s“^Mpc“^ is the 
Hubble constant. 

The subhaloes in these si mulations are identified using 
SUBFIND dSpringel et al.ll2 001 ), with merge r trees built with 
the DTrees algorithm ( Jiang et al.ll20141 . Most results in 
this paper are expected to be not sensitive to the choice of 
halo finder or tree builder. The mass, m, of a subhalo is de¬ 
fined to be its self-bound mass. The exact definition of infall 
could have some ambiguity, depending on how the bound¬ 
ary of the host halo is defined, and also because a subhalo 
may cross the boundary several times during its orbit. To 
avoid this ambiguity, we define the infall mass, rriacc, to be 
the maximum bound mass a subhalo ever had in all previ¬ 
ous snapshots. Adopting a different definition could lead to 
slightly different parameter values but does not qualitatively 
affect any of the assumptions or conclusions of our model. 
The current position of a subhalo is defined to be the current 
position of the most-bound particle it had when it was last 
resolved. In this way we follow both resolved subhaloes and 
those stripped below the subhalo mass resolution (defined 
to be 20 bound particles). 


3 THE MODEL IN A NUTSHELL 

The basic features of our model can be demonstrated in a 
simplified version as follows. Let us assume: 

(i) The unevolved subhalo mass function, i.e., the distri¬ 
bution of infall masses, rriacc, of the subhaloes accreted at 
all previous redshifts, is a power-law function, 

dN .„ 

1 n ^ rriacc 7 

d In rriacc 

with a. ~ 0.9 according to previous results fe.g. lGiocoli et aP 

l2008h . 

(ii) The unevolved spatial distribution of subhaloes, i.e., 
the spatial distribution with a given infall massQ traces the 
mass density profile, p{R), of the host halo, 

dA(i7| rriacc) 


d^R 


oc p{R), 


where d^R = dnR^dR is the volume element. 

(iii) The mass of the subhalo evolves due to tidal stripping 
after infall. The stripping is stronger at a smaller radius, 
which can be parametrized as 

- oc R^. 


This means the final mass, m, is fully determined by rriacc 
and R. For isothermal density profiles, tidal stripping pre¬ 
dicts /3 = 1. For more realistic density profiles such as NFW, 
we expect /3 ~ 1. 

Following these assumptions, the evolved distribution is 
given by 

dN{m,R) _ d A (rriacc, 7?) din rriacc 
dlnmd^i? dlnrriaccd^A dlnm 
_ d A (rriacc) d In rriacc ,, 
d In rriacc d In m 

oc m~°‘R'^p{R), (1) 

with 7 = af3 ~ 1. 

Equation [T] immediately reveals three very interesting 
features of the subhalo distribution: 

(i) The final subhalo mass function shares the same slope 
as the infall mass function. 

(ii) The final subhalo number density profile is shallower 
than the host halo density proHle (or equivalently, the un¬ 
evolved subhalo number density profile), by a factor R?. 

(iii) The subhalo mass function and spatial distribution 
are separable, meaning that the spatial distribution of dif¬ 
ferent mass subhaloes are the same except for a change in 
amplitude. 

These features agree qualitatively with existing results on 
the subhalo distribution. In Fig. [T1 the predicted spatial pro¬ 
file is compared with the simulations, focusing on its shape. 
The ratio between the subhalo spatial prohle and the host 
halo density profile indeed agrees very well with a power-law, 
with 7 = 1.3 for Aquarius halo A. 

The difference between the unevolved and evolved sub¬ 
halo spatial distribution can be understood as the result of 
a selection function, which is illustrated in Fig. [2] According 
to our first two assumptions, the spatial profiles of subhaloes 
at fixed infall mass all have the same shape, with more mas¬ 
sive subhaloes having a profile with a lower normalization. 
In the presence of mass stripping, subhaloes selected with 
the same final mass correspond to populations with differ¬ 
ent infall mass at different radii. Those found in the inner 
halo are on average more stripped, so they started with a 
larger infall mass. Hence the factor by which the profile 
is suppressed describes how the amplitude of the unevolved 
spatial profile varies with infall mass for subhaloes selected 
to have the same final mass at different radii. 


4 VERIFYING MODEL COMPONENTS 

Below we carefully test and adapt the assumptions of the 
basic model presented above using simulations. We focus on 
halo A from the Aquarius simulations as our prime example. 


4.1 The spatial PDF of accreted objects 


5 This is not to be confused with the spatial distribution of sub¬ 
haloes at the infall time. By definition, at the infall time the 
subhaloes are all located near the host halo boundary, while the 
unevolved spatial distribution refers to the distribution of the fi¬ 
nal positions of subhaloes in each infall mass bin. 


The density profile of the host halo can be interpreted as 
the probability distribution function (PDF) for the current 
position of particles accreted by the halo. As a first approx¬ 
imation, we assume the same PDF applies to each accreted 
subhalo, regardless of its mass at accretion. This PDF is 

















4 


J. Han et al. 




Figure 1. Spatial distribution of subhaloes in AqAl. The data points with Poisson errorbars show the number density profile of subhaloes 
resolved with more than 1000 particles (or 10~®M200, where M 200 is the virial mass of the host halo) at 2 = 0. The black solid line 
shows the matter density profile of the host halo. Both profiles are normalized by their values at -R200? lii® host virial radius. Right: ratio 
between subhalo number density and the DM density of the host halo, normalized to unity at -R200- The points are from the same data 
as in the left. The green solid line adopts the best-fit Einasto profile of ISpringel et al] ll2008al') for the subhalo number density. In both 
panels, the red dashed line show a power-law fit of the form (-R/-R20o)^ to PSub/pHalo inside -R200- 



Figure 2. Illustration of the spatial profile of subhaloes. Grey 
solid lines are the spatial profiles of subhaloes with a given infall 
mass, mace, with thicker lines corresponding to larger ma.ee- The 
red dashed line is the profile for subhaloes with a given final mass, 
m. For subhaloes with same m, those found in the inner halo have 
larger mace because they are more stripped. The number densities 
of subhaloes selected with different mace at each radius form the 
profile of the subhaloes with the same m. 


then given by the normalized density profile of the host halo 


p{R) 


pjR) 

M200 ’ 


( 2 ) 


where p{R) is the density profile of the host and M 200 is the 
virial mass of the host halo. Adopting an NFW profile for 
the host 


{R/R ,){1 + R/R,Y' 

where ps = 1// yi+R/R normalization, and 

Rs is the scale radius. 

The above assumption is natural if subhaloes follow sim¬ 
ilar orbits to dark matter particles. In a steady-state halo, 
the density profile is fully determined by the distribution of 
orbits of the particles, because the distribution of particles 
arou nd their orbits is fixed by the travel time at each posi¬ 
tion (|Han et al.ll201!Tl 'l. If subhaloes follow a similar distribu¬ 
tion of orbits to that of the particles, they would naturally 
form the same spatial profile as that of the dark matter 
particles. The dynamics of subhaloes differs from that of 
particles due to dynamical friction. However, for small sub¬ 
haloes for which dynamical friction is not important, one 
would expect a quite similar PDF to that of DM particles. 
Indeed it has been found that subhaloes have a distribu¬ 
tion of orbits similar to th at of the particles , with a very 
weak mass dependence (e.g. Ijiang et al.|[2015h . This means 
subhaloes are, in terms of their dynamics and spatial distri¬ 
bution, approximately indistinguishable from each other, as 
well as from dark matter particles, as long as each subhalo 
is persistently traced. Such a picture can be summarized as 
“unbiased accretion” of subhaloes relative to the accretion 
of dark matter particles. 

In Fig. [ 3 ] we test this assumption directly with the 
Aquarius simulations. We plot the density profile of the host 
halo and the number density profile of subhaloes, both nor¬ 
malized by their density at the host virial radius i? 200 - The 
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Figure 3. Left: Density profiles of subhaloes in Aquarius halo A with infall mass rriacc > 10~^M200i where M 200 is the host halo virial 
mass at 2 = 0. The profiles are normalized by their value at -R 20 O 1 the host virial radius. Different solid lines correspond to results from 
different resolution simulations of the same halo (A1 to A5), each selected with the same infall mass threshold. This threshold translates 
into different threshold in the number of particles at accretion, N^cc, for each simulation as labelled. For comparison, the dashed line 
show the density profile of the host halo, si milarly normalized by its value at i?200' The vertical short lines in the bottom mark the 
convergence radii at the five resolution levels llNavarro et aljr201ol i. Right: Ratio of the subhalo number density profile to the host halo 
mass density profile in Aquarius halo Al, normalized to unity on large scale. Subhaloes are binned according to their infall mass rriacc. 
Errorbars are shown for the lowest and highest mass bins according to the Poisson noise in the subhalo counts. The horizontal dotted 
line marks psub/pHalo = 1 for reference. 


subhaloes include all the objects that have been accreted 
by the host halo at any previous redshift, selected accord¬ 
ing to their mass at the time of infall, rriacc • The positions of 
these objects are defined by the current position of the most- 
bound particle from when the subhalo was last resolved. Ex¬ 
cept for a cuspier inner profile due to dynamical friction, the 
subhaloes closely follow the DM distribution. This is seen 
more clearly in the right panel, where subhaloes in halo Al 
are split into infall mass bins. The subhaloes closely follow 
the DM profile of the host halo, except for an excess in the 
very inner part which is more prominent for more massive 
subhaloes. The cuspier inner profile is consistent with expec¬ 
tations from dynamical friction, which causes subhaloes to 
sink toward the halo centre and is stronger for more massive 
objects. For the majority of subhaloes which are located in 
the outer halo and the low mass end, however, this effect is 
not important. 

One may question the reliability of using the most 
bound particle to extract the radial distribution for un¬ 
resolved subhaloes. The choice of using the most bound par¬ 
ticle to represent unresolve d subhaloes is in line with some 
semi-a nalytical models (e.g. ISoringel et al.|[200ll : iGao et ^ 
[20 o 3). In these models, orphan galaxies are associated with 
the most bound particles, and the resulting radial distribu¬ 
tion reproduces that of observed galaxies well. It also closely 
follows that of the DM. 

The reliability of tracking unresolved subhaloes with 
most bound particles can also be checked directly in the 
simulations. In the left panel of Fig. |3]we carry out a con¬ 
vergence study of the unevolved profile. As the resolution 
of the simulation increases from AqA5 to AqAl by ~ 5 or- 



Figure 4. The profile of subhaloes with infall mass rriacc > 
10“^M2 oo in Al, decomposed into resolved and unresolved (or¬ 
phan) populations at 2 = 0. The profiles are normalized by the 
total subhalo number density at R200- For reference, the dark 
matter density profile normalized to unity at R 200 is also shown 
as the thin solid line. 


ders of magnitude in particle numbers, more and more unre¬ 
solved objects become resolved. However, we do not observe 
any significant change in the unevolved subhalo distribu- 





































6 J. Han et al. 




^/^200 


Figure 5. Number density profile for resolved subhaloes. This is similar to Fig. [T] but only tor subhaloes that are still resolved at 2 = 0, 
with infall mass rriacc > 10“^M200. Left: subhalo profiles resolved at different resolutions in halo A, commonly normalized by the total 
subhalo number density at R 200 in Al. The black solid line is the DM density profile, normalized to unity at R200- Right: the ratio 
between the subhalo profiles and the DM profile as in the left panel. Poisson errorbars are only shown for the Al curve for clarity. The 
thick solid line is that for all accreted subhaloes in Al. 


tion. Note that even at our highest resolution, only half of 
the subhaloes in the left panel of Fig. [3] remain resolved at 
2 = 0. These resolved objects dominate at large radii, with a 
profile that still follows that of DM as seen in Fig. [I] A con¬ 
vergence study of the spatial distribution of these resolved 
objects is carried out in Fig.[Sl As the resolution of the simu¬ 
lation increases, more and more subhaloes can be resolved at 
2 = 0, and the shape of the subhalo profile approaches that 
of the DM profile down to smaller and smaller radii. The 
amplitude of the subhalo profile also converges to ~ 60% of 
that of the entire accreted population as seen in the right 
panel of Fig.[3 This suggests that about ~ 40% of accreted 
subhaloes are completely disrupted, and cannot be resolved 
no matter how much the resolution increases. This disrup¬ 
tion fraction is independent of the radial position within the 
host halo. We will discuss further support to this disrupted 
fraction when studying the stripping PDF in the Section fOl 


4.2 The unevolved subhalo mass function 

The abundance of subhaloes accreted by the host over all 
redshifts, at each infall mass mace, is known as the unevolved 
subhalo mass function (Ivan den Bosch et al.lbOOsI f. In prin¬ 
ciple, this mass function can be obtained semi-analytically 
from EPS merger trees, or analytically by combining the pro¬ 
genitor mass function predicted by EPS t heory with model s 
of a universal halo mass accretion history (lYang et al.ll201lh . 
or from EPS alone by co nsidering t he environmental depen¬ 
dence of halo formation (lLedl2004li . In Fig. [6] we show the 
unevolved subhalo mass function for halo Al, which is well 
fit by a power law 


diV 

acc 


M200 

mo 



(4) 



Figure 6. The unevolved subhalo mass function of halo AqAl. 
The red curve is the data and the solid straight line is a power-law 
fit in the form of Eq. [4] 


where mo = 10^^h~^MQ is the mass unit and A and a are 
parameters. The M 200 factor is mot ivated by the scal ing of 
the evolved subhalo mass function llGao et al.ll2004lt . For 
halo AqAl the best-fitting parameters are a = 0.96 and 
Aacc = 0.072. 


dlnm. 
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4.3 Tidal stripping of subhaloes 
4-3.1 Theoretical motivation 


Tidal stripping reduces the mass of a subhalo as it orbits 
inside the host halo. A simple model for the current mass of 
the subhalo is obtained from the mass inside its tidal radius 
rt- For spherically symmetric density profiles, the instanta¬ 
neous tidal radius can be obt ained by solving the following 
equation (see e.g.. iKiii 3 119621, for an equivalent expression) 


p(rt) 


p{R) 


2 + — 


MR) ' 

-p{R) . ’ 


(5) 


where p(rt) is the average density of the subhalo inside rt, 
R is the radial location of the subhalo, p{R) is the average 
density of the host halo inside R, and p{R) is the density 
of the host at R. The parameter Kt = Vt/Vc{R) arises due 
to centrifugal forces in a frame corotating with the subhalo 
around the host, where vt is the transverse velocity of the 
subhalo, Vc ~ GM{R)/R, and M{R) is the host halo mass 
within radius R. Once rt is known, we can calculate the in¬ 
stantaneous bound fraction p — mlm^^ of the subhalo at 
R, which we name the stripping function fJ,{R). In principle, 
u{R) depends on the density profile parameters of both the 
host and the subhalo. For isothermal density profiles, com¬ 
bined with the virial definition of masses, it is easy to show 
that n{R) only depends on the host halo size through 


p{R) 


1 R 
y /1 + Kt R200 ’ 


( 6 ) 


where R 200 is the virial radius of the host halo. Equivalent, 
but more complicated, calculations can be performed for 
NFW profiles. However for simplicity and motivated by the 
isothermal case, we model the average stripping function of 
realistic haloes with the approximation 


with a slope ~ 1, and a normalization Note the power- 
law form of p{R) is not a requirement for our model. Any 
other form can be used without affecting our conclusions. 

In reality, the stripping of subhaloes is not instanta¬ 
neous. More importantly, subhaloes are mainly stripped dur¬ 
ing their pericentric passages and once stripped the mass 
is not full y regained as the subhaloes move to the outer 
halo fe.g.. lllan et al.ll2012bl '). As a result, the current mass 
of a subhalo also depends on its past trajectory, and can be 
substantially smaller than the predicted mass inside the in¬ 
stantaneous tidal radius. So the prediction from Eq.[5]should 
be interpreted as an upper limit to the current mass. Despite 
this, the average scaling is roughly a scaled version of the 
upper limit prediction, which, as we will see shortly, can be 
well approximated by Eq. 0. Some scatter around the av¬ 
erage scaling is also expected, due to the different profiles 
and trajectories of subhaloes and an evolving host halo pro¬ 
file. The profile of subhaloes can also be modified by tidal 
heating in addition to truncation. Subhaloes can also fall 
into the host hierarchically as a subhalo of another subhalo, 
so that its current mass is also shaped by the tidal field of 
its host subhalo. All these complexities further add to the 
scatter around the average stripping function. To be more 
precise, we will use Eq. 0 to model the median stripping 
function and model the scatter around it with a log-normal 
distribution. 


For NFW density profiles, it can be shown that a sub¬ 
halo becomes completely unbound once its size is stripped to 
below 0.77rs (iHavashi et al.|[2003l f assuming the bound mass 
does not redistribute after stripping. Motivated by this, we 
expect a certain fraction of the accreted subhaloes to be 
physically disrupted at present, irrespective of the numerical 
resolution. As a result, the fraction of surviving subhaloes 
saturates as we move to higher and higher resolution, as was 
already seen in Fig. [S] 

Combining the above arguments, we model the current 
mass distribution of accreted subhaloes as a mixed distribu¬ 
tion of the following form 

dP{m\mi,cc, R) = (1 — fs)S{m)dm+ 

fsAf(\n-l7l—,\np{R),a]dlnm. ( 8 ) 

\ ^acc / 

The first term describes the physically disrupted population, 
where fa is the survival rate (1 — /s is the fraction of physi¬ 
cally disrupted subhaloes) and S{x) is the Dirac delta func¬ 
tion. The second term represents a scaled log-normal distri¬ 
bution, where J\f{x, x, a) is the probability density function 
of a normal variable of x with mean x and standard devia¬ 
tion a, and p(R) is the average stripping function. 


4-3.2 Validation of the stripping model 

In Fig. [71 we show the distribution of surviving subhaloes at 
a given radius in the different resolution simulations of halo 
A. Subhaloes are selected with the same infall mass cut in 
different simulations, in order to select the same population 
of objects. Since subhaloes can only be resolved down to a 
certain number of particles corresponding to mass m*^, the 
mass fraction sample is only complete down to 
where is the lower limit of the infall mass. As a result, 
the completeness limit in p differs across simulations of dif¬ 
ferent resolutions, and we separate each line into complete 
(solid) and incomplete (dashed) portions. Above the sam¬ 
ple completeness limit, the distributions from the different 
resolution simulations agree well. As p approaches zero, the 
cumulative distribution saturates at fa « 0.55, as expected. 
The empirical distribution can be well fit by our model dis¬ 
tribution (Eq. 0), except for the very high p end. The de¬ 
viation from a log-normal distribution at high p is expected 
because m is constrained to be below the instantaneous tidal 
limit, or at least smaller than mace by definition. We will re¬ 
visit this limit when evaluating the evolved subhalo mass 
function in Section (5) 

In the above model we have assumed the survival rate to 
be independent of the infall mass. This is examined directly 
in Fig.jS] The drop in fa at the low mass end is due to lack of 
numerical resolution, which limits the ability to resolve the 
descendants of these small objects. The fluctuation at the 
very high mass end (mace ~ 10 “®M 2 oo and above) is likely 
to be caused by statistical noise due to the small number of 
massive subhaloes available. In the dynamical range reliably 
probed by our simulation, the data are consistent with a 
constant fa- Note the survival rate is also independent of 
radius as seen in Fig. [S] These two independences further 
support the approximation that the subhaloes behave like 
indistinguishable particles in terms of their dynamics. 

The radial dependence of the distribution of surviving 
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Figure 7. Top: the cumulative distribution of for subhaloes in 
the radial range 0.5 —O.8-R200 with infall mass ruacc > 10 “^M 2 oo- 
The thick grey line is a fit to AqAl in the form of Eq. (IHll. For 
each line, the solid segment corresponds to where the sample is 
complete in /i for subhaloes with more than 100 particles at the 
present, while the dashed segment refers to the incomplete sec¬ 
tion. Bottom: the differential version of the top panel, showing 
the probability density in In fi. 


subhaloes is studied in Fig. We only include subhaloes 
with more than 1000 particles ( 10 “®M 2 oo) at infall in this 
hgure. As discussed above, the subhalo sample is only com¬ 
plete down to a mass fraction of Two complete¬ 

ness limits are shown, corresponding to a of 20 and 100 
particles respectively. 

A typical instantaneous tidal limit is shown in the left 
panel. This is calculated from Eq. © for a satellite with 
an NFW profile at accretion, with a mass of 10“‘^M2oo 
and a concentration det ermined from the average mass- 
concentration relation of ll Ludlow et al.|[2014l l. Overall, the 
tidal limit approximately delineates the upper envelope of 
the stripped masses, in agreement with expectation. A few 
data points lie above the tidal limit, because the exact tidal 
limit depends on the exact density profile of each satellite, 
which could differ from the one assumed in the calculation. 
In addition, since tidal stripping is not instantaneous, some 



Figure 8. The survival rate fs for subhaloes above a given infall 
mass, ruacc, in halo AqAl. The vertical dashed line marks where 
the infall mass corresponds to 1000 particles, below which the 
estimate of fs is limited by the incompleteness of the sample. 

subhaloes may not have been subjected to stripping long 
enough to have lost all the mass outside the tidal radius. 

Two sets of percentiles of the distribution in fj, are ex¬ 
tracted at each radius in Fig. [9l First, the median and Tier 
percentiles are extracted from all the accreted subhaloes 
that are still resolved at present. However, such percentiles 
reveal nothing about the disrupted fraction. Complementing 
these estimates, we can extract the percentiles using all the 
accreted subhaloes. This is possible because even though we 
do not know the exact /r value for unresolved or disrupted 
subhaloes, we do know their abundance at each radius and 
expect them to lie below the current completeness limit, 
above which we can still extract percentiles for the full sam¬ 
ple. In the presence of disrupted objects, the p-ih percentile 
of the surviving population would correspond to a p'-th per¬ 
centile of the full sample where p'/ 100 = (1 — fs) + /sp/100. 
This is indeed what we see in Fig. [O] Below the complete¬ 
ness limit, however, the two sets of percentiles diverge. In 
this regime, the existence of unresolved subhaloes leads to 
our hrst method overestimating the percentiles. On the other 
hand, unresolved subhaloes would cause our second method 
to underestimate the percentiles. As a result, we expect the 
true percentiles to lie somewhere in between the two esti¬ 
mates. 

Before htting the percentiles, we can first test an im¬ 
portant prediction of our model. According to our simpli¬ 
fied model (Eq. [TJ , the ratio of the subhalo count profile to 
host halo density prohle, psub/puaio, is directly shaped by 
the stripping function, fl{R), so that psub/puaio oc 
We will see later that it is also expected in a more complete 
model. To test this prediction, we scale (pSub/pHaio)^^“ ob¬ 
tained from Fig. [T] to match the amplitude of the median 
stripping function in Fig. [9] Above the completeness limit, 
the radial dependence of this ratio matches the two esti¬ 
mates of the stripping function well. Below the completeness 
limit, it is encouraging to see that the calibrated prohle ratio 
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Figure 9. The bound fraction of subhaloes in halo Al. Left: The cyan dots in the background show the location and bound fraction of 
individual subhaloes. The red solid line marks the median bound fraction of resolved subhaloes at each radius while the red dashed lines 
mark the 16th and 84th percentiles (i.e., dzlcr confidence intervals). The black lines mark the corresponding percentiles when un-resolved 
and disrupted subhaloes are also considered. The magenta solid line in the top is the expected bound fraction inside the instantaneous 
tidal radius for a subhalo inside a MW sized halo according to Eq. (0. The red dotted lines mark the limits above which the sample 
is complete in p, for subhaloes with more than 20 and 100 particles at the present time respectively. Right: Same as the left panel, 
but also showing model fits. The blue circles with errorbars are the Psub/PHalo profile in Fig. |3] tilted by l/« and scaled to match the 
normalization of the median bound fraction. The green lines are the fitted median and percentiles of the resolved subhaloes. The black 
and red lines are identical to those in the left panel. 


lies safely in between the two estimates, consistent with our 
expectation of the true stripping function. 

To obtain a parametrized stripping function, we fit a 
powerlaw to the scaled density profile ratio, (psub/pHaio)^^“. 
This fitted line is then shifted by a constant = 1.1 verti¬ 
cally which matches well the Tier percentiles. The stripping 
function is almost independent of the infall mass, as can be 
seen from Fig. [Tn]which h as a different infa ll mass selection 
to that in Fig. [9] (see also IXie &: Gad[2015l l. 

The above agreements demonstrate that our model 
adopting a power-law stripping function with a constant log¬ 
normal scatter is quite consistent with the current data. 


5 A STATISTICAL MODEL 

The above analysis reveals that our first two assumptions 
are valid and the third one is also mostly correct except 
for the scatter around the average stripping function. In the 
presence of scatter, a more complete treatment starts from 
the joint distribution of m, mace and R: 



Figure 10. Same as Fig.|9l but restricted to subhaloes with more 
than 10000 particles at infall. The green lines (model) are identical 
to those in Fig. |9] 


diV(m, mace, R) = diV(macc)/5(i?)dP(m|macc, R), (9) 

where dA’(macc), piR) and dP(m|macc,P) are given by 
Eqs. Q, ([2l) and (O respectively. Marginalizing over the 


infall mass gives the distribution of stripped subhaloes: 

/= 
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where mmin and rrimax delimit the mass range of the pro¬ 
genitors that contribute to the final population of objects of 
mass m. i? is a normalization factor arising from the integral 
over the log-normal distribution between mmin and mmax: 



erf 


( Mmin) “t“ 




\/2o 


— erf 


1 

2 


1 — erf 


/ ln(^/^max) + ao" 

V y2a 


ln(^//rmax) + aa^ 


y/2o 


( 11 ) 


( 12 ) 


with /imin = mlrrims.^ « mlM 2 oo and pmax = 

The second equality holds when m ^ M 200 so that /imin ^ 
1. A reasonable estimate of pmax is the tidal limit, which 
is roughly proportional to m(-R) as seen in Fig. making 
the second term of B a constant. For typical values of a 
and a, adopting /Xmax/M = 5 ~ 10 yields B « 0.7 ~ 0.9. 
As noted in Section ISl the log-normal distribution does 
not describe well the high /r tail of the actual distribution 
shown in Fig. 0 and so the precise value of B could differ 
from the above estimate. By matching the predicted subhalo 
abundance to simulation data, we find a universal value B ~ 
0.6 in both cluster and galaxy haloes, corresponding to an 
effective fimax/fi = 4.2. We fix B to this value hereafter. 

Setting cr = 0 and /a = 1 recovers the simple model in 
section |31 When cr 7 ^^ 0, the cr-dependent term can be un¬ 
derstood as arising from Eddington bias in selecting mass m 
subhaloes from the infall population with mace- This would 
further modify the radial dependence of the subhalo profile 
if a depends on R. We find that a constant a ~ 1 is a good 
approximation to the data. In principle, /a could also be 
radially dependent. However, because we are mostly inter¬ 
ested in the distribution of surviving subhaloes, our model 
still holds as long as the spatial PDF of surviving subhaloes, 
fap{R), follows that of the host halo. This is a good approx¬ 
imation as seen in Fig. [5] Fig. [9] indicates that /a could be 
larger at smaller radii. This could compensate for the cus- 
pier inner profile of accreted subhaloes in Fig. [3] leading to 
a profile of surviving subhaloes shadowing that of the dark 
matter. 

This distribution is still separable in m and R, leading 
to a mass-independent spatial distribution of subhaloes as 
well as a position-independent subhalo mass function. As in 
the simplified model, the evolved subhalo mass function is 
predicted to have the same slope as that of the unevolved 
subhalo mass function. We check for this directly in Fig . 1111 
where we show the ratio between the evolved and unevolved 
subhalo mass functions. For any given resolution, it appears 
that the evolved mass function has a shallower slope than 
the unevolved one at the low mass end. However, this dif¬ 
ference is strongly resolution dependent. As the resolution 
increases, the ratio becomes more and more consistent with 
a constant over the entire mass range, leading to a conclu¬ 
sion that the underlying fully resolved mass function would 
have the same slope as that of the un-evolved one. In princi¬ 
ple, a small difference can still be introduced by a weak infall 
mass dependence of fj,{R) or a deviation from the log-normal 
distribution of the stripped fraction p. 

In Fig. [T^ we calculate the predicted profile according 
to Eq. (IIOII . with parameters obtained from the previous 
sections. A good match is found between the model and the 
data, for different selections in subhalo mass. This means 



m/Maoo 

Figure 11. The ratio between un-evolved and final subhalo mass 
functions, in different resolution simulations of Aquarius halo A. 



Figure 12. The spatial distribution of subhaloes in halo Al. The 
points are data from the simulation, selected with increasing lower 
mass cuts from top to bottom. The dashed lines are model profiles 
with corresponding cuts. 


both the subhalo mass function and the subhalo spatial dis¬ 
tribution are simultaneously reproduced by the model. 


5.1 Host mass dependence 

Focusing on halo A in the Aquarius simulation, we have 
demonstrated above that a tidal stripping model coupled to 
a stochastic subhalo profile can describe well the final distri¬ 
bution of subhaloes. In this section we extend the analysis 
to a sample of haloes from the Aquarius and Phoenix simu¬ 
lations, focusing on the mass dependence of the model. The 
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Figure 13. Stacked subhalo profiles. This is analogous to Fig. 1121 
but showing the stacked distribution of subhaloes in the Aquarius 
and Phoenix simulations. The solid lines are the data and the 
dashed lines are the model predictions. The three sets of lines 
correspond to stacking subhaloes with m/M2oo above 10“^, 10“® 
and 10“® respectively from bottom to top. 


equations describing the model components apply equally 
well to these haloes. To facilitate general applications of the 
model, we list in Table [T] the model parameters extracted 
from these simulations. For the Phoenix simulation, we have 
only used the first six haloes which are close to each other 
in mass. 

The values for most of the parameters are quite similar 
between the Aquarius and Phoenix haloes. For low precision 
(e.g., at a level of ~ 10% in subhalo counts) applications, it 
is sufficient to simply take the average parameter values as 
universal and ignore any dependence on halo proHle param¬ 
eters. For applications requiring higher precision, one can 
interpolate the parameters with respect to halo mass, for 
A, and j3, as listed in the “Joint” row of Table [T] We 
leave more sophisticated calibration of these parameters in 
different haloes to future work. 

As a consistency check, the stacked subhalo radial dis¬ 
tributions in both sets of simulations are shown in Fig. 1131 
As before, the model correctly reproduces the spatial dis¬ 
tributions for subhaloes of different masses, in both sets of 
haloes. 

Analytical manipulation of Equation © may not al¬ 
ways be easy to handle. However, as Eq. (0) is given as a 
distribution function, it is easy to make Monte-Carlo realiza¬ 
tions of subhaloes in the parameter space oim, R and mace- 
Once generated, many other distributions involving these 
quantities are straightforward to evaluate. For this purpose, 
we have provided an example Python code, SuBGENOthat 
generates fast subh alo mocks according to Eg. (1^ . making 
use of the emcee (iForeman-Mackev et al.l l2013l l Markov- 


Chain-Monte-Carlo sampler. We will also be taking a Monte- 
Carlo approach in some of the following applications. 


6 DISCUSSION & APPLICATIONS 

The model has important implications for modelling galaxy 
formation and the statistics of subhaloes in simulations 
themselves. Once calibrated, the model can be applied to 
predict many population properties of subhaloes and galax¬ 
ies. We give some examples below. 


6.1 Do galaxies trace subhaloes? 


The connection between the unevolved and evolved spa¬ 
tial profiles has clear implications for understanding the 
connection between subhaloes and galaxies. Taking halo- 
occupation distribution models as an example, an important 
ingredient in such models is the spatial proHle of member 
galaxies that are used to populate haloes with galaxies. It 
is common practice to populate haloes with galaxies follow¬ 
ing the density profile of the host halo. This choice is mo¬ 
tivated either due to simplicity by assuming galaxies trace 
dark matter or by the galaxy dis tribution in semi-an alytical 
models or SPH simulations (e.g.. iBerlind et ^l2003h . If one 
assumes galaxies trace subhaloes, however, it appears that 
the spatial distribution of subhaloes from numerical simu¬ 
lations conflict with this model because the spatial profile 
of subhaloes appears flatter in the inner halo than that of 
dark matter. The same discrepancy exists when comparing 
the observed radial distribution of galaxies with the sim¬ 
ulated distribution of subhaloes, a phenomenon termed as 
the anti-bias of subhaloes. This discrepancy is erased if one 
recognizes the differences in the two profiles is an effect of 
sample selection, or an improper comparison between galaxy 
and subhalo samples. The flatter profile in simulations is a 
result of selecting in bound mass. In observations, galax¬ 
ies are more likely to be selected according to stellar mass, 
which corresponds more closely to selecting subhaloes by in¬ 
fall mass. In this case, planting galaxies following the host 
halo profile is the corre ct choice. Such a selection effect has 
alread y been noted b y Nagai fc Kravts^ (l2005| i: Kravtsov! 
(l2010ll . I Conroy et al.l ll2006[ ) and IChen et al.l (l2006l ~l showed 
that populating subhaloes according to their infall proper¬ 
ties provides a better match to the observed spatial distri¬ 
bution or clustering of galaxies, although they only studied 
surviving subhaloes (i.e., no orphan galaxies). 

Note the anti-bias is not purely a result of failing to 
model orphan galaxies. The sample selection bias also oper¬ 
ates when only surviving subhaloes are considered. 


6.2 The “universal” subhalo mass function 

Subhalo mass functions, both unevolved and evolved, are 
shown in Fig. 1141 The evolved subhalo mass function can be 
obtained by marginalizing over the spatial part of Eq. m. 
giving 


dN M200 r 'm 

— -= A 200 - - 

a mm mo [mo 
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Table 1. Parameters of the full model, extracted from the Aquarius (Milky Way) and Phoenix (Galaxy cluster) simulations. For each 
simulation set, the parameters are estimated from the stacked distribution of the level 2 haloes. The “Joint” row provides interpolated 
parameter values intended for application to arbitrary halo masses, mo = is simply the adopted mass unit. The last row 

lists the equations in which each parameter is defined. 



M2oo/(/i-iMo) 

Aacc 

a 


P 

a 

fs 

Aquarius 

(1.0 ±0.3) X I 012 

0.089 

0.95 

0.42 

1.4 

1.1 

0.54 

Phoenix 

(6.7 ± 1.0) X lOi'^ 

0.080 

0.95 

0.34 

1.0 

1.1 

0.56 

Joint 


0.1(M20o/mo)“°-‘’^ 

0.95 

0.5(M20o/'mo)“°'°® 

1.7(M2oo/mo)-°'°'‘ 

1.1 

0.55 

Reference 


Eq. l|4} 

Eq. ijUl 

Eq. 10 

Eq. 0 

Eq. 0 

Eq. 0 


with 

(14) 

Consistent with the good agreement seen in Fig. 1131 the 
predicted subhalo mass functions match the data very well. 
Snbstituting the parameter values and the numerical halo 
density profiles into Equation m, we obtain approximately 
a common normalization parameter ^200 — 0.008 for both 
sets of simulations, meaning the model automatically repro¬ 
duces the roughly un iversal subhalo m ass function known in 
previous studies (e.g. ICao et aI1l2004bll . For comparison, the 
dotted line in the lower regions show a joint power-law fit to 
both simulations, which is almost identical to the individual 
predictions. A combined fit to the unevolved mass function 
is also shown in the upper region by a dotted line, which 
unsurprisingly takes the average parameter values listed in 
Table [T] 

The model predicts an abundance ratio A 2 oo/Aacc — 
0.1. It is interesting to note that this is the same as the typ¬ 
ical mass fraction locked in subhaloes. The connection be¬ 
tween the two would be straightforward if we assume that we 
can extrapolate the unevolved subhalo mass function down 
to smooth accretion, ruacc —>■ 0, or if the majority of the 
halo mass is accreted from subhaloes. In both cases, the host 
halo mass can be found by adding up contributions from all 
the progenitors, M 200 = / m-accdVacc, so the subhalo mass 
fraction is simply / mdN/M 2 oo = A 2 oo/Aacc. We have ex¬ 
plicitly checked that in Aquarius Al the resolved accretion 
from subhaloes adds up to 90% of the host halo mass. 

The above subhalo mass function is defined to be the 
abundance of subhaloes inside R 2 oo- A better understanding 
of the universality of this mass function can be obtained if 
we generalize it to be defined inside a different radius, as 


dN{< R) ^ 

d In m mo mo 


(15) 


The normalization A{R) can be predicted from our model 
through 


A{R) = AaccB/sC' 


P{r)n{R) “dV 

M(< R) 


(16) 


In absence of stripping, i.e., for ft(i?) = 1, the subhalo dis¬ 
tribution simply follows the host density profile. Then the 
normalization factor is independent of R, and also largely in¬ 
dependent of the host halo mass as determined by Aacc/s. In¬ 


side a real halo, however, stripping is important, and the spe¬ 
cific abundance of subhaloes would depend on both R and 
the host halo mass. This is shown in Fig. 1151 For R > R 200 , 
tidal stripping is not important. In addition, the majority 
of the subhalo population inside R is comprised of those 
close to R. This results in the abundance A{R) being largely 
determined by the infall population, leading to an approxi¬ 
mately universal specific abundance A across different halo 
masses. For R <C R 200 , however, both a mass and radius de¬ 
pendence is introduced by tidal stripping, because the sub¬ 
halo profile no longer follows the host profile. The universal¬ 
ity of the subhalo mass function thus can be understood as a 
consequence of the fact that subhaloes trace the density field 
unbiasedly on large scales {R > i? 2 oo)- It is good to see that 
two other definitions of the virial radius, i?vir which is de¬ 
fined to enclose an averag e density accord i ng to the spherical 
collap se prediction fe.s.. I^e et al.l[l996l : iBrvan fc NormanI 
Il998ll , and i? 200 b defined to enclose an average density of 200 
times the mean matter density, both lie beyond i? 200 - Thus 
the subhalo mass function defined inside Rvir and i? 200 b are 
also approximately universal. 

Note the specific abundance plotted in Fig. [TS] is cal¬ 
culated from the simulation data rather than from Equa¬ 
tion (0. This is because Equation GH) is not accurate for 
R » R 200 where the integration over the log-normal distri¬ 
bution in Equation (El) breaks down (so that B is no longer 
a constant) due to the constraint mlm^cc < 1. However, 
our Monte-Carlo sampler SubGen can handle this bound¬ 
ary condition easily and predicts the right subhalo abun¬ 
dance and spatial distribution even outside iZ 200 - 


6.3 Dark matter annihilation emission from 
subhaloes 

The dark matter annihilation emission from subhaloes in 
massive haloes is a ve ry promising target to search for dar k 
matter particles (e.g., iGao et aI1l2012al : IPinzke et al.ll20li] ~l . 
With our current model for the distribution of subhaloes, we 
can proceed to calculate the annihilation emission integrated 
over the subhalo population. This is achieved by modelling 
each subhalo with a truncated NFW profile. The NFW pa¬ 
rameters are fixed according to the mass and concentration 
at infall, while the truncation radius can be obtained from 
the current mass of the subhalo. The relevant equations to 
obtain the annihilation emission from such a truncated pro¬ 
file are listed in Appendix If we combine our model with 
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Figure 14. The subhalo mass functions. The shaded regions 
above and below are the unevolved and evolved subhalo mass 
functions respectively, with the dashed lines marking the av¬ 
erage mass function and the regions spanning the 1-cr scatter. 
The red and green parts are the results from the Aquarius and 
Phoenix simulations respectively. The upper two solid black lines 
are power-law fits (Eq. ((Jl)) to the unevolved mass function with 
parameters listed in the first two columns of Table ^ The lower 
two black solid lines are subsequent model predictions of the 
evolved mass functions. The dotted lines are joint fits to the two 
simulations bearing the average parameter values. 



-^200 


Figure 15. The specific abundance of subhaloes, defined as the 
normalization of the subhalo mass function in Equation lISli. 
The two lines show the average abundances in the Aquarius and 
Phoenix haloes respectively. The three vertical lines correspond 
to three common definitions of the virial radius from left to right 
respectively: R 2 Q 0 , Rvir ^ 200 b (see text for definitions). 


Figure 16. The mass-concentration relations (at red shift 0) used 
in this work llMaccio et al.ll2008l : iLudlow et 51l2014ll . The green 
shaded region shows the = 0.3 scatter which we apply to 

both models. 

a mass-concentration relation for the subhaloes at infall, we 
can obtain both the spatial distribution of subhaloes and the 
truncated density profile of each subhalo. This enables us to 
evaluate the annihilation emission of each object as well as 
their distribution inside the host, down to any mass limit of 
interest. Such a calculation is easy to do with our subhalo 
sampler, SubGen, and can reproduce the key features of the 
spatial and mass dependences of the annihilation emission. 

We assume the concentrations follow a log-normal dis¬ 
tribution at fixed mass, with a scatter uiogc = 0.13. For 
the median mass-concentration relation, we consider two 
2 = 0 models from the li terature. One is a power-law 
fit from iMaccid et al.l (|2008l ) for virialized hal o es, an d the 
other is a physical model from iLudlow et al.l ll20l4) that 
calc ulates concentration from the mass accretion history (see 
alsolMa ccid et aT1|2008l : IZhao et al.ll2009l : IPrada et al.ll201^ : 
ICorrea et al.ll2015l . for some similar models). In both models 
we adopt the same cosmology as that of our simulations. As 
shown in Fig. 1161 the two models are consistent with each 
other for massive haloes resolvable by contemporary cosmo¬ 
logical simula tions. For sma l l obje cts (M 200 < IO^/i'^Mq), 
however, the ILudlow et al.l l|2ni4l ) model predicts a much 
lower concentration than that extrapolated from a power- 
law relation. 

Fig. \T7\ shows the luminosity profile of annihilation 
emission from subhaloes in a cluster sized halo. For com¬ 
parison, we also show two models gene ralized from fitting 
the Aquarius and Phoenix simulations dPinzke et al.l 120111 : 
iGao et ^ l2012al ). in normalized coordinates x = R/R 200 
and Zsub = Lauh/LsubiR < R 20 o)- Although both profiles 
are obtaine d from fitt i ng eith er Aquarius or Phoenix simu¬ 
lations, the lGao et al.l (l2012al ) model is fitted to the surface 
brightness profile (i.e., the projected luminosi ty profile) of 
subhalo emission, while the lPinzke et al.l |20l3) model is ini¬ 
tially given in the form of the cumulative luminosity profile 
Lsub(< R)- The two different ways of fitting and the differ- 
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ent parametrizations are likely to lead to disagreement where 
the profiles are poorly constrained by the simulation. This is 
evident in Fig.[T7]for R < 0.1i?2oo. Outside this region, it is 
encouraging to see the shape of our predicted luminosity pro¬ 
file is in good agreement with both models. Inside this region 
our model predicts a cuspier profile than the two existing 
models which are essentially unconstr ained by the simula - 
tions at such small radii. Adopting the lLudlow et af] (l2014l l 
profile results in less subhalo emission, while the shape of 
the luminosity profile from subhaloes is barely affected. 

The amplitude of the subhalo emission is typically spec¬ 
ified in terms of a boost factor, defined as the total annihila¬ 
tion luminosity from subhaloes normalized by that from the 
smooth component of the host halo, both evaluated inside 
the host virial radius. This is ex amined direct l y in th e right 
panel of Fig. 1171 Adopting the iMaccio et al.l (l201Cll l mass- 
concentration relation, the boost factor scales with muni ap¬ 
proximately as a power-law function, with a slope consistent 
with the value —0.226 estim ated from the resolved subhaloes 
in the Aquarius simulation (ISpringel et al.|[2008bl l. Down to 
an Earth mass, a nominal free-streaming mass scale for cold 
dark matter haloes, the boost factors rise to a few hundred 
and a few thousand for galaxy and cluster haloes respec- 
tiv ely. These v alues a re slightly higher than th ose estimated 
in I Gao et ahl (l2012al '): ISpringel et al.l ll2008bl ~l by extrapo- 
lating the Aquariu s and Phoenix simulations. When the 
iLudlow et al.l ll2014ll mass-concentration relation is adopted, 
however, the 6(miim) function is no longer a power-law, and 
is significantly reduced at low mum, reflecting the greatly 
reduced concentration of haloes at the low mass end in this 
model. Down to an Earth mass, the boost fac tor is reduced 
by a fa ctor of 50 in both haloes when using the lLudlow et al" 


{2014) relation compared with that using the lMaccio et al 
(201C) relation. 

The lowered boost factor, in addition to the cuspier 
emission profile from subhaloes, makes the total luminos¬ 
ity p rofile insid e a halo less extended than that expected 
from iGao et ^ (l2012al l ; iPinzke et al.l (120111 1 . This implies 
that constraints on the dark matter annihilation cross- 
section in clusters based on previous boost-factor estimates 
dHuang et al ] |2012l : [Han et al.l|2012al . e.g.,) could be relaxed. 
We provide some fitting formulae for the subhalo emission 
in Appendix IB] 

Our approach differs from some previous est i mates 
fe.g.. IStrigari etlih 20071 : Anderhalden fc Diemandl l2013l : 
ISanchez-Gonde &: PradcJ 2014ll in that we start from the in¬ 
fall mass to infer the density profile, rather than from the 
current mass which has been affected by tidal stripping. The 
concentration of subhaloes plays a vital role in this esti¬ 
mate, with lower concentrations leading to lower boost fac¬ 
tors. We acknowledge several limitations of our current es¬ 
timate. First, the mass-concentration relation at infall time, 
instead of that at z = 0, should be applied to the infall 
mass. This causes the concentrations to be over-estimated 
when the 2 = 0 relation is used, leading to an over-estimate 
in the boost factor. For example, lowering the concentra¬ 
tion parameters by 0.2 dex (roughly corresponding to the 
mass-concentration relation at 2 ~ 2) leads to a reduction 
in the boost factor by a factor of 3. The correct concen¬ 
tration distribution can be found by looking at the redshift 
distribution of the progenitors either in simulation or from 
EPS theory. We restrain ourselves from calibrating such rela¬ 


tions in this work. We notice that iBartels &: Andol ll2015h has 
recently combined analytical models of the unevolved sub¬ 
halo mass, the accretion-redshift distribution and a redshift 
dependent mas s-concentration with the average mass strip¬ 
ping rate from I Jiang fc van den Boschl (|2014| ~I for an evalu¬ 
ation of the boost factor. Secondly, the infall mass function 
is extrapolated to low mass with a power-law form, while 
in principle it could differ from a power law and could be 
calculated theoretically with the EPS theory. Thirdly, the 
stripping function is also assumed to be independent on in¬ 
fall mass down to the lower mass limit of subhaloes. While 
this is a very good approximation for the subhaloes resolved 
in our simulations, deviations from a powerlaw form could 
become important once the infall mass range becomes much 
larger. A simple estimate utilising Eq. ((SJ suggests a steeper 
stripping function for low mass subhaloes, which could re¬ 
duce both the boost factor and the subhalo emission in the 
inner halo. Despite these limitations, the current model still 
improves over previous work and can be extended using the 
current framework. In its current form, our predicted boost 
factors should be taken as upper limits. 


6.4 The lensing mass profile 

The mass of subhaloes as a function of projected cluster¬ 
centric radius at a fixed stellar mass can be measured with 
weak lensing (e.g. El et al.1 12OI3I. |2014 I2OI5I: ISifon et al.1 
I2OI5II . Because stellar mass is most directly related to the 
infall mass of its host subhalo, such measurements essen¬ 
tially constrain the mass of subhaloes selected at fixed infall 
mass. When stacking subhaloes, disrupted ones (assuming 
their galaxies persist) contribute no signal, thus they only 
act to dilute the average signal from surviving subhaloes. 
In the presence of disrupted subhaloes, the measured signal 
AEobs = /sASreai where AE is the difference between the 
cumulative and differential surface density profile of dark 
matter halo. For subhaloes with (truncated) NFW profiles, 
the lensing signal AE is proportional to m (when the other 
parameters are fixed). Failing to model the disrupted sub¬ 
haloes would lead to under-estimating the subhalo mass by 
a factor fa. Note this fa is not simply the fraction of or¬ 
phan galaxies in semi-analytical models, because the latter 
is not a physical quantity but depends on the resolution of 
the simulation used by the model. 

Once the disrupted fraction has been corrected for, the 
measured subhalo mass as a function of projected halo¬ 
centric distance can be obtained. Because the surviving mass 
is not a single value at a fixed infall mass, the measured mass 
is some average of the underlying mass distribution which in 
general lies in between the mean and median of the distri¬ 
bution dMandelbaum et al.ll200^ ~l. To interpret this “lensing 
average” and to correct it to the true median or mean masses 
requires knowledge of the underlying mass distribution. For 
subhalo lensing at a given projected distance, the relevant 
distribution is 

P(m|macc, 7?p) = [ P{fi\R)p{R)dl. (17) 
J 1.0.S 

Given this distribution, the mean and median subhalo mass 
can be evaluated analytically or, more conveniently, with the 
Monte-Carlo sampler SubGen. The generated Monte-Carlo 
samples can also be used to evaluate the systematic bias in 
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Figure 17. Left: The luminosity density from subhaloes in a cluster sized halo (M 200 = 6.7 x 10^^h~^M(T)) . The green solid line and 
the red dashed line are the predicted annihilation emission from our model adopting the lLudlow et ahl ll2014h and iMaccio et ^ ll2008l^ 
mass-concentrations respectively. For c omparison, the blue solid line show the annihilation emission from the smooth density distribution 
of the host halo. Two model profiles llGao et al.l 120123 : IPinzke et al.] 1201 if) that extrapolate the annihilation emission from resolved 
subhaloes in simulations are also shown. All the luminosities are normalized by the total luminosity from the smooth density distribution 
of host halo. Right: the boost factor contributed by subhaloes above a given mass limit miim- The red and green lines are for the galaxy 
and clus ter sized halo respe ctively. For each halo, the dashed line show the result when the a power-law avera ge mass-concentrat ion 
relation llMaccio et alj[2010l) is adopted, and the solid line show that when a physical mass-concentration model llLudlow et al.ll2014h is 
adopte d. The red and green arrows on the vertical axis mark the extrapolated boost factors for the two haloes according to | Gao et alJ 
|20123) down to earth mass (lO^^M©). The short grey line shows a power-law scaling with b oc llSpringel et al.ll2008bl) . 


the lensing measurement relative to the real median mass, 
by simulat i ng th e lensing averaging process as was done in 
iHan et al] (120141 ') for real observations. 

In Fig. 1181 we compare our predictions with a recent 
measu rement of subhalo mass in galaxy clusters by lLi et al.1 
1 I 2 OI 5 II . To populate subhaloes with galaxies (i.e., converting 
mace to m*), we adopt th e stell ar mass-infall mass relation 
determined in IWang et al.l (l2013l l for satellite galaxies with a 
scatter criogm. = 0.19 at Hxed infall mass. We have corrected 
for the different definit ions between our infall mass and that 
in IWang et al.1 (|2013| ). The subhaloes are further selected 
with a stellar mass threshold to compare with the obser¬ 
vations. Our result i s very similar to that from the mock 
catalogue in iLi et ^ (l2015l ) created from a semi-analytical 
galaxy formation model, but requires much less effort to ob¬ 
tain. Accounting for the disrupted subhaloes increases the 
measured mass by a factor of ~ 2. After this correction, 
the measurements start to show a significant tension with 
the predicted mean and median. However a full investiga¬ 
tion which would have to consider many issues is beyond the 
scope of this paper. For example, the observational selection 
function is more complicated than w e have assumed an d in¬ 
volves selection in host halo mass llSifon et al.ll2015h and 
redshift. Systematic uncertainties in stellar mass estimates 
may also introduce bias in the mass ratio as well as com¬ 
plicating the selection function. Contamination from neigh¬ 
bouring massive groups is likely to cause an over-estimate 
at large radii. On the other hand, the value of /s in the real 
universe may differ from the value used here. Our /a is es¬ 
timated from dark matter only simulations. The value of /a 


in the real universe may be different due to baryons, which 
make subhaloes more resistant to tidal disruption. By ap¬ 
plying the /a correction we have also assumed that galaxies 
are not disrupted together with their host subhaloes, which 
may not be the case in the real universe. High resolution 
hydrodynamical simulations are required to address these 
uncertainties. 


7 SUMMARY & CONCLUSIONS 

We have developed a model that unihes the distribution of 
subhaloes in mass, m, position, R, and infall mass, mace. 
The model fully specihes the joint distribution of these three 
quantities in an analytical form (i.e. Equation 0: 

dY(m,macc,R) = diV(macc)p(i?)dP(m|macc, R), 

where dY(macc) describes the infall mass distribution, p{R) 
is the spatial probability distribution of dark matter parti¬ 
cles inside the host halo, and dP(m|macc,P) describes the 
hnal mass distribution of subhaloes of a given infall mass 
at a given radius. The specihe forms of the relevant terms 
in the joint distribution are given by Equations @, © and 
([5|l , with parameter values applicable to different host haloes 
masses listed in Table [T] A Monte-Carlo sampler, SubGen, 
is also provided that easily generates subhalo samples inside 
any host halo following the above distribution. Once a sub¬ 
halo sample is generated, any population statistics involving 
these variables can be easily obtained. 





























































16 J. Han et al. 


a median radial dependence well approximated by a power 
law. 



Figure 18. The projected profile of subhalo mass to stellar mass 
ratio in galaxy clusters. The dashed and solid lines represent the 
mean and median mass of survived subhaloes with stellar mass 
M* > in a cluster with M 200 = The 

shaded region is bounded by the 15th and 75th percentiles of the 
sub to stellar mass ratio at each radiu s. The ci r cles w ith error- 
bars are the original measurements from iLi et al.l ||2015^ while the 
triangles are original results multiplied by l//s to account for the 
disrupted subhaloes. 


The support for this model can be summarized as fol¬ 
lows: 

• Using high resolution ACDM cosmological simulations 
of b oth a galaxy and a c luster sized hal o from the Aquar - 
ius llSpringel et al.l [2008al l and Phoenix llGao et al.ll20129 l 
projects, we have carefully verified that the shape of the un¬ 
evolved spatial distribution (i.e., the radial profile at fixed 
mace) follows the density profile of the host halo, a phe¬ 
nomenon we summarize as unbiased accretion of subhaloes. 
This holds for both surviving subhaloes and unresolved or 
disrupted subhaloes as traced by their most-bound parti¬ 
cles. Dynamical friction leads to a deviation of the unevolved 
spatial distribution from that of the host halo density pro¬ 
file only in the very inner region and is important only for 
subhaloes with very large macc/Af 2 oo- 

• The amplitude of the unevolved spatial distribution, 
as described by the unevolved subhalo mass function, 
dA’/d In mace, follows a power law in each individual halo. 

• The joint distribution is then obtained following Bayes 
theorem, by further specifying the connection between m 
and mace with the conditional distribution P{m\ms.cc, R)- 
This connection is shaped by tidal stripping, with subhaloes 
in the inner halo being more heavily stripped on average. 
Through a convergence study, we find that about 45% of 
subhaloes are physically disrupted (i.e., stripped to m = 0 
regardless of numerical resolution). Because the spatial dis¬ 
tribution is independent of infall mass, the same disruption 
fraction applies to all infall masses and at all radii. For the 
surviving subhaloes, we find P(m|macc,^Z) can be approx¬ 
imated by a log-normal distribution at each radius, with 


Marginalizing (i.e., integrating) the joint distribution 
over any variable, one obtains the joint distribution of the 
remaining ones. For example, marginalizing over the infall 
mass, the model simultaneously reproduces the universal fi¬ 
nal mass function and the universal spatial distribution of 
subhaloes of a given final mass. In particular, the model 
predicts that: 

• The final mass function follows a power-law form with 
the same slope as the infall mass function. 

• The spatial distribution of subhaloes at fixed m, which 
we call the evolved spatial distribution, is flatter than the 
density profile of the host halo. The ratio between the two is 
determined by the amount of tidal stripping at each radius. 
This explains the so called “anti-bias” between the galaxy 
distribution and the subhalo distribution as purely a selec¬ 
tion effect. 

• The shape of the evolved distribution is also indepen¬ 
dent of m. The scale-free nature (i.e., power-law form) of 
the infall mass function and the mass-independence of the 
unevolved spatial profile are the keys to such independence. 

The parameters of our model ingredients have been cal¬ 
ibrated with simulations and we find only very modest vari¬ 
ation with simulated halo mass. This enables the model to 
be safely interpolated to other halo masses. The calibrated 
model can be applied to a wide range of problems. We give 
several such examples, including the universality of the sub¬ 
halo mass function, the dark matter annihilation emission 
from subhaloes, and lensing measurements of subhalo mass. 

We demonstrate that the universality of the subhalo 
mass function exists because subhaloes trace the density 
field at large radii where tidal stripping is irrelevant. Further 
inside this radius, the mass function is lower in more massive 
haloes. Using the framework to calculate the dark matter 
annihilation emission of subhaloes, we demonstrate that the 
adopted mass-concentration relation for subhaloes is crucial 
in such calculations. Extrapolated down to an Earth mass, 
the commonly adopted powerlaw mass-concentration model 
overpredicts the total subhalo emission by a factor of 50 
compared with the results when adopting a more physical 
mass-concentration relation. The model can also be easily 
adapted to compare with, as well as to calibrate, gravita¬ 
tional lensing measurements of the subhalo mass. The ex¬ 
istence of a physically disrupted subhalo population could 
potentially lead to a correction to the lensing measurement 
by a factor of ~ 2, amplifying the tension be tween a recent 
subhalo lensing measurement llLi et al JImi) and theoreti¬ 
cal predictions. 

The model can be extended to higher redshift and fur¬ 
ther calibrated in other cosmologies. A dependence on host 
halo concentration may also be introduced as additional 
model parameters. The aspect of the model that is of the 
most interest and least known is how subhaloes are stripped. 
This is described in the model by the average stripping func¬ 
tion and its scatter. The unevolved subhalo mass function, 
on the other hand, can be fully predicted from EPS theory. 
In addition, EPS calculations are also capable of providing 
the distribution of accretion redshifts, which can be com¬ 
bined with a redshift-dependent mass-concentration relation 
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to provide accurate density profile parameters for subhaloes. 
This could for example improve the predictions for the sub¬ 
halo annihilation emission. 
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APPENDIX A: ANNIHILATION EMISSION 
FROM A TRUNCATED NEW HALO 


The mass and concentration parameters of an NFW pro¬ 
file can be easily converted to (ps, ra), the scale density and 
scale radius of the halo. If the profile is truncated at n, the 
truncated mass, m = m{rt), can be related to Vt through 


m = rriB 


ln(l + xt) 


Xt 

1 + Xt 


n = - 


1-b 


Wo (—e“P+™'Gt)/»7Js)) 


(Al) 

(A2) 


where nis = 47rpsrs, Xt = n/rg. Wo{x) is the principle 
branch of the Lambert W function. Subject to a factor that 
is determined by the particle properties of DM0 the anni¬ 
hilation emission of a truncated NFW profile is 


-L(77T., TUacc, Cacc) — 


)= rp\r)d^ 

Jo 


1-(1+^)-^ 

Vs 


47r 2 3 

3 " 


(A3) 


This luminosity depends very weakly on the truncation ra¬ 
dius when Vt > Xa- 


® Note our results in Section 16.31 and Appendix El are always 
expressed in terms of the normalized luminosity L/Lhaiot which 
are independent of the particle physics factor. 


APPENDIX B: FITTING FORMULAS TO THE 
SUBHALO ANNIHILATION EMISSION 


The annihilation emission from subhaloes is usually factor¬ 
ized as 


dLsub(-R) T T clZ/sub(^) 

d^R d^x 


(Bl) 


where Lhost is the total emission from the smooth density 
held of the host halo, x = R/R 200 is the normalized radius, 
Lauh{x) = Lauh{x)/Lauh{x = 1) is the normalized luminosity 
with Laub( 2 ;) being the total emission from subhaloes inside 
X, and h is the boost factor so that La,A{x = 1) = feLhoat- 
An analytical function that hts our luminosity proHles 
reasonably well is 


dLsub ( 3 :) 
d®* 


= 0.1(a;-b0.15)' 


(B2) 


This prohle is not sensitive to the mass of the host halo. 
It can be further projected to obtain the surface brightness 
prohle for observational applications. The boost factors from 
subhaloes above an Earth mass can be htted with 


b = 4.6 


M 200 


(B3) 


^/i-1Mq 

when adopting the lMaccio et al] (l2010l l mass-concentration 
and 


b = 0.08 


M 200 


(B4) 


>-iM0, 

when the lLudlow et all (l2014h relation is used. 

For reference, the lumin osity prohle and bo ost factors 
from iGao et ^ (|2012a|) a nd Pinzke et al.l (l201lll are listed 
below. For the lGao et al.l ll2012ah model, 


di)sub(:r) 

d^a; 

b 


4.53(1 -b , 

1.6 X 1O-®(M2oo/M0)°-®® 


(B5) 

(B6) 


Equation l[B5l) i s obtaine d by de-projecting Eq.(2) in 
iGao et ^ J2012al in For the I Pinzke et al.l (120111 ) ht, 


dZsub(a;) _ ai[l-b 02 ln(a;)]a:“i''“^+“2 


d^a; 


6 = 0.017 


4:TTX^ 

M 200 


lO-6h-iM0 


(B7) 

(B8) 


with ai = 0.95 and 02 = —0.27. Equation (iBTll is ob- 
tained from differen tiating the original cumulative profile in 
IPinzke et al.l (|2Qlll) and an earth mass of has 

been adopted in Equation (|B8|) . 


^ There is a factor of 2.5 difference in normalization between 
Eq. (|B5|l and the original de-projected version. This is because 
we normalize the profile by t he tot al luminosity inside the 3-D 
R < R 200 , while iGao et ahl l|2Q12al ') normalize it by that inside 
the projected radius Rp < -R200- 


























































